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Abstract. Calculations of the final merger stage of binaiy black hole evolution can only be carried 
out using full scale numerical relativity simulations. We review the status of these calculations, 
highlighting recent progress and current challenges. 


INTRODUCTION 

The final coalescence of a binary black hole (BBH) is driven by gravitational wave 
emission and proceeds in 3 stages [1]. In the quasi-adiabatic inspiral phase, the BHs 
have a large enough separation that they can be treated analytically as point particles 
in the post-Newtonian limit. This relatively slow inspiral is followed by a dynamical 
merger phase, during which the BHs plunge toward one another and merge, forming 
a single BH. This highly distorted remnant then emits gravitational radiation in quasi- 
normal modes as it rings down to form a quiescent rotating BH. 

Calculating the gravitational radiation signatures from BBH coalescence is essential 
for understanding and interpreting the data expected from both ground-based and space- 
based detectors [2]. Both the inspiral and ringdown stages can be handled analytically. 
During the inspiral, the waveform is a “chirp,” which is a sinusoid increasing in both 
frequency and amplitude as the orbital separation shrinks. Since the detectors can typi- 
cally observe thousands of orbits during inspiral, these waveforms can be used as tem- 
plates for data analysis algorithms based on matched filtering to yield information on the 
masses and spins of the components. The final ringdown of the distorted remnant will 
produce a burst of gravitational radiation. The waveforms from this stage are damped 
sinusoids that can be calculated analytically using BH perturbation theory, and can be 
used in identifying the properties of the final Kerr hole. 

However, the merger phase, starting at the end of the inspiral and proceeding through 
the plunge and dynamical merger, is governed by strong, nonlinear gravitational fields. 
The resulting burst waveforms can only be calculated through numerical relativity sim- 
ulations of the full Einstein equations in 3 spatial dimensions (3-D) plus time. Gravita- 
tional wave observations of BH mergers can provide an outstanding probe of a regime 
expected to be rich in strong-field spacetime dynamics, including black hole spin flips 
and couplings. 

Numerical modeling of BH mergers a challenging enterprise, due to the deep mathe- 
matical nature of general relativity, the need to solve a fairly large number of nonlinear, 
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coupled partial differential equations, and the wide range of physical scales involved. 
While we still do not have definitive waveforms starting in the late inspiral stage and 
proceeding through ringdown, significant progress has been made over the past few 
years. In this article, we provide a brief review of the status of numerical relativity BH 
simulations, highlighting recent successes and current challenges. While this inevitably 
involves some excursions into the mathematical side of general relativity, this will be 
done for purposes of illustration rather than rigor. Further details and additional refer- 
ences for interested readers can be found in the recent review by Baumgarte and Shapiro 

[3]. 


EINSTEIN’S EQUATIONS IN 3+1 FORM 


Provided that any accompanying accretion disks are dynamically negligible, BBHs 
coalescing through the emission of gravitational radiation are solutions to the vacuum 
(i.e., source-free) Einstein equations of general relativity. As such, their timescales are 
proportional to the total mass M and, after time scaling, all other properties of the 
dynamics and waveforms depend only on ratios involving the masses and spins of the 
components. This has the pleasant consequence that calculations of the gravitational 
waveforms from all three stages of BBH coalescence can be easily scaled to apply to 
any masses and spins [4]. 

The Einstein equations in vacuum are [5] 

(4) ^v"5 (4) ^v = 0, (1) 

where we use spacetime indices ju, v = 0, 1,2,3. Here, is the 4-D or spacetime 

Ricci tensor, R = ( 4 )R |iV g /iV , and we assume an implied summation over pairs of 
repeated up and down indices. Since all second rank tensors used here are symmetric 
(e.g., = ^R Vfl ), Eqs. (1) comprise 10 coupled, nonlinear partial differential 

equations for the evolution of the spacetime metric g^ v which are second order in both 
time and space derivatives. 

To evolve BBH mergers using numerical relativity, we must first write these equations 
in a form suitable for solution on a computer. We will accomplish this using the so- 
called 3+1 spacetime split, in which 4-D spacetime is envisioned as being sliced into a 
collection of 3-D spacelike hypersurfaces, threaded by a congruence of timelike curves 
along which time t is measured [5]. 

The coordinate freedom of general relativity provides 4 freely specifiable coordinate 
conditions or gauge choices. In the 3 + 1 framework, these are the lapse function a, 
which measures the lapse of proper time adt between neighboring slices, and the 3 
components of the shift vector (3‘, which governs the movement of spatial coordinates 
during evolution from slice to slice. The spacetime metric takes the form 


8[iv — 


-a 2 + ft/3* ft \ 

Pj % ) ' 


( 2 ) 
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where y i; - is the spatial metric of the 3D slice, /3,- = y 1 */ and spatial indices i,j= 1,2,3. 

We use units in which c = G = 1. The choices made for a and j5‘ are critical for 
successful BH evolutions, as discussed below. 

The Einstein equations (1) split into 2 sets. Four equations set conditions on the 
spacelike slices: the Hamiltonian constraint 

R + K 2 -K u K ij = 0 (3) 

and the momentum constraints 

D J K J i -D i K = 0. (4) 

Here, R = y^/J^. is the Ricci scalar and Dj is the spatial covariant derivative compatible 
with the 3-metric y, ; - in a spacelike slice. K- is the extrinsic curvature of the slice, and 
K = K\ = y' 7 X; * s * ts h'ace. The constraint equations (3) and (4) also set conditions 
on the initial data, as discussed in the next section. The remaining 6 equations give the 
evolution of the initial data in time. Taking y,- • and K to be the basic variables, we can 
write them as a set of equations that are first order in time. In the standard Amowit- 
Deser- Misner (ADM) [6, 3] formalism, these become the evolution equations for the 
3 -metric y^-, 

d,Y,j = -2 ctKy+Dfij+Dfi, (5) 

and for the extrinsic curvature K t j, 

S,l fy = -D,Dja + a(R lt - 2 K^j + KKy) + P k D t K y + K„pfi k + K kj Dfi k . (6) 

Note that the terms involving and R contain both first and second spatial derivatives 
of Yij- 


INITIAL DATA FOR BBH MERGERS 

Numerical simulations of BBH mergers should begin with initial conditions characteriz- 
ing the late stages of the inspiral, before the plunge begins. Since the binaries reach this 
state due to the emission of gravitational radiation, realistic initial conditions must take 
account of this radiation as well as the positions and momenta of the BHs. Techniques 
are available for producing BBH initial data that satisfies the Einstein equations, and 
current efforts are focused on insuring that the data represents the astrophysical state of 
the binary accurately. 

Setting up initial data for a numerical relativity simulation requires solving the 4 
constraint equations (3) and (4) on a spacelike slice for y f • and [7]. Most efforts to 
solve these equations begin with a conformal decomposition of the metric and extrinsic 
curvature. Specifically, the physical 3-metric y-y is written 

Yij = W%, (7) 

Calculating Gravitational Wave Signatures from Binary Black Hole Mergers July 4, 2003 3 



where iff is known as the conformal factor and f t j is the conformal metric. The extrinsic 
curvature is split into its trace K and a trace-free part Ay, so that 

K U =A U + frjK. ( 8 ) 

Various approaches then differ by the way in which A - is decomposed [8]. 

In all of these techniques, the initial data variables separate into 2 sets: those variables 
that are freely specifiable (such as the conformal metric y the trace K, and some parts 
of Ay), and those that are constrained (typically, the conformal factor iff and other parts 
of Ay). Once the freely specifiable quantities are chosen, they are inserted into the 
conformally transformed initial value equations; these nonlinear elliptic equations are 
then solved for the constrained parts of the data. Combining the free and constrained 
pieces then produces a valid set (yy,/sTy) of initial data that obeys the Einstein initial 
value equations Eqs. (3) and (4). 

In most of the work on the initial value problem to date, the freely specifiable data have 
been chosen to decouple and otherwise simplify the solution of the constraint equations. 
For example, we can choose the initial slice to be conformally flat (jf.j = the metric of 
flat space) and have K — 0 (known as maximal slicing). With a particular decomposition 
of Ay, this not only decouples the constraint equations but also allows analytic solutions 
of the momentum constraints for a BH with both linear momentum and spin [9]. Since 
the momentum constraints become linear equations for the constrained parts of Ay in 
this approach, we can superpose these solutions to obtain data for any number of BHs 
with linear and angular momentum. The nonlinear Hamiltonian constraint must still be 
solved numerically. This type of data, pioneered by Bowen and York [9] has formed the 
starting point for various numerical simulations of BHS, including grazing collisions 
[10] confirm!! and merger calculations beginning near the final plunge [11, 12]. 

Much current work centers on producing BBH initial data sets that are astrophysically 
realistic. In general, this involves choosing spatial metrics that are not conformally 
flat. Some recent work has attempted to incorporate the results of post-Newtonian 
expansions near the end of the inspiral phase into the freely specifiable data (e.g., [13]) 
while other efforts have used a non-flat conformal 3-metric based on superposing 2 
BHs in Kerr-Schild coordinates (e.g. [14]). A special challenge arises because, in the 
conformal decomposition process, the physical metric y- and the extrinsic curvature 
K t j that emerge out of the solution process can be changed significantly from the 
freely specifiable quantities that are chosen at the start. The conformal thin-sandwich 
decomposition [15], which includes information on the evolution of the metric away 
from the initial slice, may alleviate such problems. It is also a promising approach for 
constructing BBH initial data on quasi-circular orbits near the plunge [3]. 

All of these techniques can introduce unphysical gravitational waves into the initial 
data. For example, a recent study constructed BBH initial data sets using similar choices 
for the freely specifiable data but different conformal decompositions [8]. The amount 
of gravitational radiation present differed among the resulting initial data sets by a few 
percent of the total mass, an amount that is comparable with the total gravitational 
radiation expected from the entire merger process. 
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EVOLVING BBH MERGERS 


The BBH initial data must be evolved for ~ 3.5 orbits or more near the end of the 
inspiral, through the plunge and merger, and into the ringdown stage [16]. For a BBH 
system near plunge with separation a ~ 6 M, the orbital period P ~ 90 M. We therefore 
need to be able to evolve the system for ~ 1000M or longer, to produce astrophysically 
relevant waveforms useful for detection. 

Over the past decade, considerable effort has been expended by the worldwide nu- 
merical relativity community towards achieving this goal. While significant progress 
has been made, a number of difficult challenges remain. In this section we’ll present a 
brief overview of BBH evolutions, highlighting both current successes and outstanding 
problems. 


Choice of Formalism 

The ADM evolution equations (5) and (6) mark the starting point for most BH 
simulations. However, a direct numerical implementation of these equations in 3-D 
quickly develops exponentially growing unstable modes that cause the code to crash. 
In particular, a BBH evolution using the standard ADM equations lasts only ~ 13M 
[1 1], a small fraction of an orbital period near the start of plunge. While the exact cause 
of these instabilities is not yet understood, we do know that they are properties of the 
mathematical formulation of the Einstein equations, and not of the numerical methods 
employed [17]. 

Shibata and Nakamura [18] and Baumgarte and Shapiro [19] modified the original 
ADM evolution equations (5) and (6) using a conformal-traceless decomposition. This 
so-called BSSN formalism is obtained by writing the conformal factor in the form 
y/ = e^, so that 

fij = e ~ 4<t> Yij, (9) 

with the choice f = det(^y) = 1. The traceless part of K- is scaled according to 

A IJ = e- 4 *A lj . (10) 

The conformal connection functions 

P EE f>% = ~3jf, (11) 

are introduced as new independent variables. Here, the Tj k are the connection coeffi- 
cients (or Christoffel symbols) associated with fjj, and the second equality in (1 1) relies 
on the condition f = 1. The ADM evolution equations then take the form 

40 = ~ s a ^+/3 l 4'0 + 54^ l > (12) 

d t Yij = -2 aAij + P k d k Yij + Y ik dfi k + Y kj dfi k - (13) 

d t K = —y ij DjD i a + a(A i:j A ij + \K 2 ) + )3 i d i K, (14) 
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d t Aij 

d t P = 


= (-(D,.D y a) rF + a(Rffj) + a{KA tj - 2 A a A l j) 

+P k d k A ij -{-A ik djP k +A kj dfi k - \ A {j d k P k , 

-lA^djCC + 2a (P jk A kj - l? j djK + 6A'^) 
+&djP - ridjp + fr^.jSl + \fdjdfii + fj^djl 3\ 


(16) 


where the superscript TF in (15) indicates the trace-free part (e.g. = 7? ( y - y- ; 7?/3) 

and the definition (11) serves as a new constraint equation. Using f‘ and its first spatial 
derivatives, the Ricci tensor is then written so that the only second spatial derivatives 

of that remain are in the Laplacian f lm d l d m f i j', see [19, 3] for details. 

The BSSN equations (12) - (16) constitute a set of 17 evolution equations. In a 
typical simulation, all of these equations are updated, and the conditions y = 1 and 
trace (A,- -) = A*- = 0 are then imposed. One or more of the new constraints (11) may also 
be imposed. 

Numerical implementations of the original BSSN formalism given here exhibit much 
better stability properties than those using the standard ADM equations [20]. Further 
improvements and variations have been introduced, mostly based on adding some of 
the constraint equations to the evolution equations, or enforcing some of the constraints 
during evolution [21, 22] (see also [23]). With these developments, single rotating BHs 
can now be evolved stably for several thousand M [22] and BBHs for > 100M [24]. 

A somewhat different approach to developing improved formalisms centers on writing 
the Einstein equations in hyperbolic form. Most current efforts along these lines result 
in large sets of equations that are fully first order in both time and space derivatives 
and have many adjustable parameters [25, 26]. The stability of these systems is strongly 
dependent on the way these parameters are chosen [17, 25], Currently, evolutions of a 
single BH are possible for up to several thousand M for optimal parameter choices [27]. 
Although these systems are mathematically elegant, the large numbers of equations and 
free parameters involved can make them unattractive from the point of view of numerical 
implementation. However, they do have the advantage that their characteristic structure 
can be analyzed readily and used to impose natural boundary conditions both at the outer 
edges of a numerical grid and on inner excision boundaries of BHs [25]. 


Representing a BH on a Numerical Grid 


A BH contains a physical singularity, at which the curvature is infinite. Numerical 
relativity simulations involving BHs thus face special challenges in representing these 
objects on a computational grid. Several techniques have been developed to accomplish 
this task. 

The spatial metric for a Schwarzschild BH in isotropic coordinates is conformally 
flat, so that y ; • = y/ 4 /^-. Taking K- = 0, the momentum constraints are trivially satisfied 

and the Hamiltonian constraint becomes A flat yr = 0, where A flat is the Laplacian in flat 
space. With the asymptotically flat boundary condition yr — » 1 as r — > the conformal 

factor takes the expected form yr = 1 +M/2r, where M is the Schwarzschild mass and 
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the event horizon is located at r H = M/2. Since the Hamiltonian constraint is a linear 
equation in this case, we can easily add solutions to get initial data for multiple BHs. 
Thus the conformal factor for 2 BHs is yr = 1 + My/2r^ + M 2 /2r 2 , where the BH of 
massMj is located at withr^ = (jc — x^ + Cy — yj) 2 + (z — Zj) 2 , and similarly 

for M 2 [28]. 

For a single BH, the coordinate transformation r 1 = M 2 /4r is an isometry that maps 
every point inside r n — M/ 2 into a point outside r H . The geometry inside r H is thus the 
same as that outside r H , allowing us to think of the BH as consisting of 2 asymptotically 
flat regions or sheets connected by a throat at r = r H [5]. (need to clarify where the 
physical singularity went) For 2 BHs, we can have a single sheet with 2 throats, each 
of which connects to a separately asymptotically flat region, giving a 3-sheeted topology. 
Alternatively, we can add additional throats, corresponding to mirror images of the 
companion BH, inside each of the original throats. This conformal imaging approach 
restores the isometry and yields a 2-sheeted topology in which the throats connect 
identical asymptotically flat sheets [29]. To implement this in numerical relativity, we 
can use the isometry conditions on the throats as boundary conditions, thereby removing 
the singularities inside from the grid cite ??. 

A similar procedure can be carried out following the Bowen- York [9] approach to 
yield data for BHs with spins and linear momenta. As discussed above, we assume 
conformal flatness and take K = 0. Then A,- ; - takes a specific form that solves the 
momentum constraints; since these are linear equations in this case, the solutions can 
be added together. The Hamiltonian constraint is a nonlinear equation that must be 
solved numerically. To remove the BH singularities from the grid, isometry boundary 
conditions have to be applied at the throats, which are rather complicated surfaces [30]. 

The puncture method eliminates the need to apply isometry conditions on the throats 
[31]. We consider the initial slice to have points, called punctures, removed at and 
r 2 . The conformal factor is written as the sum iff = u+y/', where if/' = 1 + / 2r j + 

J$ 2 /2r 2 , and the constants -> Mj and Jt 2 -4 M 2 as the separation between the BHs 
approaches infinity. The Hamiltonian constraint becomes a nonlinear equation for the 
function u. In this approach, the singular terms are absorbed into the analytic term yr' 
and the remaining function u is regular everywhere, even at the punctures r l and r 2 . The 
puncture technique has also been extended to post-Newtonian data [13]. 

Puncture data can be evolved by choosing the location of the punctures so that they 
do not coincide with a point on the computational grid. With suitably chosen lapse 
and shift, the punctures do not move through the grid (although their positions in 
spacetime, governed by the evolution of the metric, do change), and the spacetime in 
their immediate vicinity does not evolve [24]. The puncture method has been used to 
evolve single and multiple BH systems, including a grazing collision [10] and the final 
plunge [11, 12]. 

The other method that has been developed for handling BHs on a grid is known as 
excision. This approach relies on the fact that no physical information can progagate 
from inside the BH event horizon to influence the spacetime outside. Thus, it is not 
necessary to evolve the Einstein equations on the spacetime inside the horizon, and the 
interior of the BH can be excised from the numerical grid [32] 

The event horizon is the boundary between those spacetime events which emit light 
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rays reaching infinity, and those which do not. Since the event horizon is a global quanti- 
tity, it is necessary to know the entire spacetime before the location of the event horizon 
can be found. However, this is generally not possible during a numerical simulation, in 
which the spacetime is constructed by evolving from one slice to the next. In practice, 
a related concept known as an apparent horizon is used. The location of the apparent 
horizon can be calculated on each spacelike slice using y ( - and K i j. Since the apparent 
horizon is always located inside the event horizon [33], it is safe to excise zones within 
the apparent horizon. In the case of Schwarzschild and Kerr BHs, the locations of the 
apparent and event horizons coincide. 

In a typical implementation of excision, several buffer zones separate the apparent 
horizon from the excised region inside. Several different shapes have been used to define 
the excision boundary, including a cube and a sphere. An appropriate discretization of 
the Einstein equations must be used for the grid points on the excision boundary. And, 
when an excised BH is moved, points that were previously within the excised region 
become points on the computational grid. Special extrapolation techniques are used to 
populate these points with physical data [34]. 

Excision has been used to evolve both static and distorted single BHs [22, 35], single 
BHs that move across a grid [36, 34], and a grazing collision of 2 BHs [37]. 


Choice of Lapse and Shift 

A key ingredient of successful numerical relativity simulations is the choice of suit- 
able coordinate conditions. Within the 3 + 1 framework, these are specified by the lapse 
function a, which controls the slicing of spacetime, and the shift vector /3 ', which allows 
the spatial coordinates to move [38]. 

The simplest choice is geodesic slicing, in which a = 1 and j3‘ = 0. Here, the 
coordinates coincide with freely-falling observers that travel on geodesics normal to 
each spacelike slice. However, since gravity is attractive, these observers fall towards 
the BH and will generally hit the singularity in a finite time. This cause the code to crash 
before any significant dynamical evolution of the system can be accomplished. 

Singularity avoiding coordinate conditions have been developed to prevent the space- 
like slices from intersecting BH singularities. In this approach, conditions are imposed 
on the lapse function so that the resulting sices wrap up around (but do not intersect) the 
singularity. For example, the maximal slicing condition K = 0 yields an elliptic equation 
for a that must be solved on each slice. Since the solution of elliptic equations tends 
to be computationally expensive, hyperbolic, parabolic, and algebraic conditions for a 
have been developed which also produce singularity avoiding slices. The most popular 
is the “1 + log” slicing, in which a — 1 + lny [39]. 

Singularity avoiding slices with j3 ' = 0 have enabled evolutions of BHs up to t ~ 30 — 
40 M, including a BBH grazing collision [10]. However, as the evolution proceeds the 
slices become increasingly distorted, leading to large gradients in the metric. Eventually 
the resulting grid stretching produces numerical errors that cause the code to crash. 

A nonzero shift vector can be used to keep grid points from falling towards a BH. 
Elliptic, hyperbolic, and parabolic conditions have been developed for /3 1 that counter the 
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grid stretching properties of singularity avoiding slices. Using the hyperbolic T— driver 
condition for f5‘ along with 1 + log slicing has extended the running time of single BH 
evolutions to several thousand M [22], and BBH evolutions to several hundred M [24]. 


Adaptive Mesh Refinement 

Numerical relativity simulations of BBH mergers must typically solve ~ 17 or more 
evolution equations on a 3-D spatial grid. The need to provide adequate resolution of the 
multiple physical scales in these models implies that some form of adaptivity is required 
to carry out these calculations, even with today’s high performance computers. 

For example, the gravitational waves produced typically have wavelengths ~ 10 — 100 
times the scale of the BHs themselves. A sufficiently large grid is necessary so that these 
signals can propagate into the wave zone and be extracted at large distances from the 
source. The grid must also be large enough so that the boundary conditions at the outer 
edges of the grid can be applied in an effectively asymptotic region. To accomplish these 
objectives, higher resolution can be used in the central interaction region, with coarser 
resolution in the wave zone. 

In general, the BHs in a binary will have different masses. For intermediate mass and 
supermassive BBHs, the component masses can easily differ by factors of ~ 10 — 100 or 
more, pointing to the need for different spatial resolution in the vicinity of each BH. In 
the early stages of the merger, the BHs have a relatively wide separation. Ideally, each 
BH would be surrounded by an appropriately high resolution grid, with lower resolution 
in the orbital region, and even coarser resolution in the wave zone. The grid resolution 
would then change adaptively as the merger takes place to accomodate the evolving 
physical situation. 

The use of mesh refinement in 3-D numerical relativity simulations of BHs is still in 
its early stages. Fixed mesh refinement has been applied to single BH evolution [40] and 
a short part of a BBH evolution [41]. Adaptive mesh refinement also was used in the 
propagation of a source-free gravitational wave across a grid [42]. 

In general, gravitational waves produced in a BBH merger will need to cross one or 
more mesh refinement boundaries as they propagate from the higher resolution of the in- 
teraction region into the lower resolution of the wave zone. Interpolation conditions must 
be set on the numerical data at the mesh refinement boundaries to allow the waves to 
cross them smoothly. Linear interpolation conditions, typically used in hydrodynamics 
codes, result in significant spurious reflections for gravitational waves crossing refine- 
ment boundaries. The use of quadratic interpolation conditions at the mesh refinement 
boundaries dramatically reduces the spurious reflected waves [43]. 


Additional Challenges 

There are several other areas which are considered critical for the success of numerical 
relativity BBH simulations. One key area is the role of the constraint equations, and in 
particular how best to incorporate these relationships into the solution of the Einstein 
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equations [44, 45]. The role of the outer boundary conditions (that is, as r ®o) and 
their relationship to preserving the constraint equations during an evolution is also under 
investigation [46]. 

While most numerical relativity simulations to date have been carried out using fairly 
standard finite difference techniques, some groups are applying spectral methods in their 
codes (e.g. [47]). Since BBH simulations typically require high performance computing 
resources, the resulting codes must be implemented in a parallel computing environment 
with attention to speed and the efficient use of memory. 


Towards Astrophysical BBH Mergers 

While the use of numerical relativity to model BBH coalescence has proved to be 
very challenging, significant progress has been made over the past decade. Although 
important issues remain, currently available techniques are sufficient to begin the study 
of astrophysically relevant BBH merger waveforms. An important step in this direction 
was taken by combining full numerical simulations of the merger of two BHs with 
perturbative techniques to evolve the late-time state of the system. This so-called Lazarus 
approach yielded the first model of a gravitational waveform from an astrophysically 
plausible merger [1 1]. 

This method begins with initial data for a BBH near the plunge. The idea is to use a 
full 3-D numerical relativity simulation to evolve this system up to the point at which 
the BHs are close enough to be treated as a single, highly distorted BH. The evolution of 
this remnant (that is, the ringdown stage) will then be carried out using techniques from 
BH perturbation theory. In this approach, the large scale simulation is reserved for the 
strongly nonlinear part of the evolution, and simpler numerical techniques are applied to 
handle the ringdown and extract the gravitational waveform. 

The Lazarus approach was first used to evolve the final plunge of 2 identical 
Schwarzschild BHs. The result of this calculation was that ~ 3% systemSs total energy 
and ~ 12% of its total angular momentum was radiated away in the form of gravitational 
waves, leaving behind a rotating Kerr BH with rotation parameter a/M ~ 0.72 [11]. 
More recently, these calculations have been extended to include the plunge of BHs with 
spin [12]. 

So far, the full numerical relativity evolution used in these calculations lasts only 
~ 13 M, as the original ADM formalism is employed used. Since the dynamics of 
realistic BBHs mergers is expected to take significantly longer than this, these results 
must be considered preliminary. Nevertheless, they represent an important advance 
and we can expect to see the Lazarus approach used in longer-lived, more realistic 
simulations in the future. 
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